Negative Differential Friction Predicted in 2D Ferroelectric In2Se3 Commensurate Contacts

Abstract At the macroscopic scale, the friction force (f) is found to increase with the normal load (N), according to the classic law of Da Vinci–Amontons, namely, f = µN, with a positive definite friction coefficient (μ). Here, first‐principles calculations are employed to predict that, the static force f, measured by the corrugation in the sliding potential energy barrier, is lowered upon increasing the normal load applied on one layer of the recently discovered ferroelectric In2Se3 over another commensurate layer of In2Se3. That is, a negative differential friction coefficient μ can be realized, which thus simultaneously breaking the classic Da Vinci–Amontons law. Such a striking and counterintuitive observation can be rationalized by the delicate interplay of the interfacial van der Waals repulsive interactions and the electrostatic energy reduction due to the enhancement of the intralayer Se—In ionic bonding via charge redistribution under load. The present findings are expected to play an instrumental role in design of high‐performance solid lubricants and mechanical‐electronic nanodevices.

In contrast to liquid lubrications, which may fail in some extreme conditions, [27,28] structural superlubricity has emerged as a new promising remedy for the reduction of friction, [29,30] which generally originates from the effective cancelation of lateral forces with incommensurate rigid crystalline contacts [31,32] due to the lattice mismatch and thus significantly reduces the energy barriers of motion. [33][34][35] Such an intriguing phenomenon was originally predicted three decades ago in nanoscale homogeneous graphitic contacts [36] and later termed superlubricity, [37] where the friction exhibits peak and superlow feature in commensurate and incommensurate contacts, respectively. However, only recent experiments have unambiguously observed the proposed superlubricity in various 2D material incommensurate van der Waals (vdW) junctions, including homogeneous contacts of graphite nanoflake-graphite, [38] heterojunctions of graphenehexagonal boron nitride (hBN), [39] and graphene-MoS 2 . [35] More intriguingly, adhesion-dependent negative friction coefficients of chemically modified graphite was reported recently in experiment. [40,41] Note that, here the negative friction coefficient could be more rigorously expressed as negative differential friction (NDF), coined by analogy with another two very popular concepts of negative differential resistance [42,43] and negative differential capacitance. [44] Moreover, NDF was also theoretically probed in superlubric graphene-hBN heterojunctions, [45,46] graphene-graphene, [47] and ferroelectric materials. [48] However, in all the previously reported systems including graphite-graphite homo-junction, the superlubricity and NDF were sustained either by the condition of incommensurate contacts, [36,37] which may be blocked to commensurate configurations during the sliding, [37] or by the reduction of the potential corrugation by the vdW interactions in the attractive regime of the interfacial separation above its equilibrium. [46,48] Very recently, nonmonotonic interfacial friction under load is predicted in 2D crystals, such as MoS 2 , hexagonal boron nitride, and graphene bilayers. [49] 2D ferroelectric material, as a unique member of 2D material family, may offer new opportunities in superlubricity due to electrically tunable interlayer dipole-dipole coupling which is absent in other types of sliding contacts. Here, we employ first-principles calculations to predict that negative friction coefficient can be realized in 2D ferroelectric In 2 Se 3 commensurate contacts which break the incommensurate contact condition and the classic Da Vinci-Amontons law. The friction force f, measured as the potential corrugation or sliding barrier (E bar ) [33][34][35]50] as rationalized by the Prandtl-Tomlinson (PT) model, is lowered upon increasing normal load N. It is widely accepted that the interlayer potentials or sliding corrugation of 2D junctions are essentially determined by two factors, vdW and electrostatic interactions. [51][52][53] Therefore, as long as the delicate interplay of the interfacial vdW force and the electrostatic energy dependent on the dipole-dipole alignments at the ferroelectric 2D contact can be tuned to appropriate regime, negative is destined to emerge, and the condition of incommensurate contact becomes dispensable. We validated this strategy in the investigation of the sliding friction between two quintuple layers (QL) of 2D ferroelectric In 2 Se 3 [54][55][56][57] commensurate contacts with appropriate interlayer dipole-dipole alignments. Here, the vdW interactions enhance the potential energy corrugation (PEC) due to the repulsive regime of the interfacial separation subject to the dipole-alignment dependent interfacial interactions under load, in contrast to previous findings [46,48] that vdW interactions predominantly reduce the PEC in the attractive regime of the interfacial separation above its equilibrium. However, the electrostatic energies reduce the PEC by enhancing the intralayer Se-In ionic bonding due to charge redistribution under load, and the delicate interplay between the vdW interactions and electrostatic energy leads to NDF.

Properties of 1QL-In 2 Se 3 and 2QL-In 2 Se 3 Systems
First, the calculated lattice constant of In 2 Se 3 is 4.11 Å and the optimized thickness of In 2 Se 3 is about 6.8 Å, [54] in close agreement with previous calculations [54] and experimental results. [55,56] As shown in Figure 1a,b, the significantly different interlayer spacing between the Se layer and the two In layers and the in-plane centrosymmetry breaking of the 1QL In 2 Se 3 result in spontaneous out-of-plane (≈0.1 eÅ) and in-plane (≈2.4 eÅ) electric polarizations, [54] respectively.
Then, we identify the preferred high-symmetry stacking modes of the homogeneous contact. As shown in Figure 1, for 2QL In 2 Se 3 commensurate contacts, in the out-of-plane direc-  Figure S1 in the Supporting Information, respectively. As also summarized in the Table S1 (Supporting Information), the binding energies of the most stable three stacking configurations are highly correlated with their interlayer distances (d), namely, the stronger the binding of the 2QL contact is, the smaller the d becomes, as manifested by the optimized d of 2.86, 2.98, and 3.02 Å, for O-O(AP), I-O(P), and I-I(AP), respectively. Such a trend is also observed for the three low-lying counterparts.

Interlayer Sliding Energy Profile of 2QL-In 2 Se 3 Homojunctions Along [110]
Now, starting from the optimized most stable stacking configuration of O-O(AP), we analyze the energy profile of the interlayer sliding along two typical commensurate pathways, i.e., [110] and [100] directions, where the bottom Se atoms of the top QL (Se BT ) are displaced straightly from A site, over B and C, to the second-nearest neighboring A, and directly displaced from an A site to the nearest neighboring A site, as respectively shown in Figure 1c. To mimic the experimental investigations of friction as performed on other 2D materials, [35] when sliding the top QL In 2 Se 3 relative to the bottom QL, the X, Y, and Z coordinates of the bottom-two layers of atoms in the bottom QL In 2 Se 3 are fully fixed, and the energy profile of the commensurate contact is optimized by displacing the top QL In 2 Se 3 every 0.3 Å along the proposed directions with fixing X and Y coordinates of the toptwo layers atoms. [34] That is, in the constrained sliding pathways,  Figure 2a-c, respectively. Two distinct features can be observed: i) In all the three contacts, during the whole sliding period, two local maxima of energy are identified, with the highest energy maximum (E max ) occurred when Se BT is slid to around the A site; ii) With the external N increasing, the E max can be slightly pass over the A site, and the lowest energy minimum (E min ) may be shifted from B site to C, as found in the case of I-I(AP) with N = 2.2 GPa. Correspondingly, we present the activation energy barrier (E bar ) of the interlayer sliding subject to different loads in Figure 2d-f, for O-O(AP), I-O(P), and I-I(AP), respectively. Here, we define the E bar as the energy difference between the E max and E min . Dividing the E bar by the sliding distance, [64] we can estimate the friction force f during the sliding. Very interestingly, the three types of commensurate contacts possess very contrasting E bar (f) curves as a function of load N, as shown in Figure 2d-. Specifically, for O-O(AP), the E bar (f) decreases with N increasing, namely, the differential friction coefficient = ∆f/∆N <0 can be identified in a fairly large load regime up to 2.2 GPa applicable in recent experiments on 2D material contacts. [12,45] Here, we emphasize that the NDF in the present load regime is also confirmed by further calculations performed by using HSE+D3 and optPBE-vdW schemes, as detailed in Figure S2 in the Supporting Information. However, > 0 is observed when N is beyond 2.2 GPa. Note that, here the predicted overall variation of friction with normal load is estimated to be about several nN according to the definition; [34] importantly, the variation of the friction within 0-2.2 GPa reaches up to about 6.4%. Therefore, we expect that the NDF in the present work can be probed by the elaborate experiments. [38] Moreover, we have also calculated the differential friction coefficients ( = ∆f /∆N) for all the three representative systems. As shown in Figure 3a-c, it is very clear that for the case of O-O(AP), the friction coefficient is indeed negative in the load regime of 0 < N ≤ 2.2 GPa, with a maximum absolute value of 0.015 at 2.2 GPa; for the case of I-O(P), is negative (with relatively small absolute values) in the load regime of N ≤1.10 GPa; In contrast, for I-I(AP), is positive and increases up to the maximum absolute value of about 0.04 at 2.2 GPa, beyond which the calculated begins to decrease with the load, corresponding to the -to--like structural distortion. These findings provide strong evidence that the electric dipole can serve as an important degree of freedom to modulate the differential friction coefficient of In 2 Se 3 interlayer commensurate sliding, i.e., from < 0, through ≈ 0, to > 0.

Origin of the Negative Differential Friction Coefficient μ
To rationalize the above intriguing findings, we decompose the calculated E bar (f) into two terms which are contributed purely by vdW interactions (E bar−vdW ) and electrostatics interactions (E bar−ES ), [51][52][53] as shown in Figure 2g  I-I(AP), respectively. Here, the E bar−vdW is defined as the energy difference of the vdW correction energies between E max and E min , and attribute the left portion of the E bar to electrostatics interactions, respectively. As seen from Figure 2g Figure S3 in the Supporting Information show that the Ewald energy and exchange correlation potential dependent on the charge transfers and charge redistributions under load dominate the E bar reduction and negative friction. For I-I(AP), in the whole load regime of ≈0-3.3 GPa, the E bar−vdW exhibits monotonous decrease which nevertheless is insufficient to compensate the monotonous increase of the E bar−ES , leading to a reversed sign of the , i.e., > 0.

Role of Interlayer vdW Interactions Verses Intralayer Electrostatic Interactions
Here we emphasize that though NDF arises from the decreased sliding potential corrugation with increased normal load is also observed in other systems, [45][46][47][48] however, such NDF is dominated by van der Waals interactions in the attractive regime of the interfacial separation above its equilibrium. [46,48] However, in our findings, the vdW interactions increase, rather than decrease the potential corrugation with the load increasing, very probably due to the fact that the interlayer distances are reduced to the sharp repulsive regime of the interfacial vdW interactions, as supported by the smallest interfacial separations in the OO(AP) and detailed analysis in Figure S4 in the Supporting Information. These findings demonstrate the critical role of the electrostatic interactions in dominating the NDF of the ferroelectric 2D material In 2 Se 3 vdW commensurate contact. Moreover, as detailed in Figure S5 in the Supporting Information, the increase of the friction under load beyond 2.75 GPa in O-O(AP) is accompanied with a structure distortion of the top QL In 2 Se 3 in the E max state, and an to phase transition is observed beyond 3.3 GPa. Such a transition results in energy increases, due to the higher energy of the phase than , by about 0.01 eV for the 1QL In 2 Se 3 54 under the present pressure regime, as also detailed in Figure S6 in the Supporting Information. For I-O(P), such a structural distortion and phase transition begins from relatively small critical load, around 1.65 GPa. However, no significant phase transition is observed in both the E max state of I-I(AP) contact and E min states of all the three systems.
To illustrate the underlying mechanism of the specific alignment of ferroelectric polarizations in leading to distinct frictions of the present systems, we first analyze the energy band structure of 1QL In 2 Se 3 in Figure S7 in the Supporting Information, which exhibits semiconducting characteristic with an energy bandgap of 0.78 eV. Due to the out-of-plane polarization, there are positive and negative charges localized respectively on the top and bottom sides of a free standing 1QL In 2 Se 3 as defined in Figure 1a. Such a feature is also supported by the centrosymmetry breaking nature (along the z direction) of the electronic charge density of the valance band, which is mainly localized on the top surface, i.e., the starting side of the electric dipole, see Figure S8 in the Supporting Information. Correspondingly, the interfacial couplings of the three 2QL In 2 Se 3 commensurate contacts can be schematically simplified as models of two contacting surfaces carrying positive and/or negative polarization charges, as shown in the left panels in The contrast friction behaviors of the three systems can be readily understood from the distinct responses of these polarization charges at the E max configurations under load. Specifically, for O-O(AP), due to the electrostatic repulsive interactions of the negative charge themselves and the experienced electric field force applied by the electric dipole of the opposite QL, part of the surface polarization negative charges are respectively repelled to the nearby Se atom sites, as demonstrated by the charge difference analysis in Figure 4d. Such a charge accumulation (in red) on Se atoms lowers the electrostatics energies of the two individual QL components and thus reduces the E bar−ES via enhancing the Se-In ionic binding. Importantly, such a charge transfer is significantly enhanced subject to the load N increase. For examples, the E bar−vdW increases from 0.129 to 0.159 eV, however the E bar−ES decrease from −0.008 to −0.046 eV when N increases

Interlayer Sliding Energy Profile of 2QL-In 2 Se 3 Homojunctions Along [100]
Now, we briefly report that when the top QL  Figure S10 in the Supporting Information. These findings suggest that the in-terlayer sliding prefers to proceed along the present "guide rail," i.e., [110] direction, probably due to the significantly anisotropic features of the potential surface, which facilitates in realizing superlubricity. [45]

Role of External Electric Field in Tuning the Friction
To the end, taking the most stable contact O-O(AP) as a typical example, we briefly highlight that the new strategy established here in modulation of the friction by ferroelectric polarization can be further convincingly demonstrated by applying external electric field which can modify the out-of-plane polarization and/or by changing the in-plane polarization arrangement of the two QLs as well. As shown in Figure 5a, for the case of O-O(AP), under positive electric fields, along the [110] direction, the E bar (f) can be further lowered and the negative friction coefficient feature becomes more significant in the relatively low load regime (N < 2.2 GPa), which can be also rationalized by the interfacial polarization charges changes under the external electric field, as detailed in Figure S11 (Supporting Information). Moreover, as shown in Figure 5b, the sliding E bar (f) in O-O(AP) can be significantly lowered via reversing the in-plane dipole alignment to parallel arrangement, i.e., O-O(P), however leading to positive in the whole present load regime. These findings suggest that In 2 Se 3 with controllable friction is highly desirable in design of functional nanodevices wherein tunable friction is demanded.

Conclusion
In conclusion, the present study predicts that negative sliding friction coefficient can be realized in electrically polarized 2D material commensurate contacts. By selecting appropriate dipole-dipole alignments between the individual components of the contact, the raise of the sliding energy barrier contributed by vdW interactions can be thoroughly compensated by the electrostatic energy reduction, which is a spontaneous polarization charge redistribution process under appropriate load regime. As a consequence, negative friction coefficient can be obtained in the commensurate contact, and incommensurate condition is broken. Also importantly, the new strategy established here in reversion of the sliding friction coefficient by electric-dipole is also demonstrated by applying out-of-plane electric field and/or by changing the in-plane ferroelectric polarization alignments. Due to the fundamental importance, the present central findings are also expected to play an instrumental role in developing functional lubrications, nanosensors, and related mechanical-electronic nanodevices based on ferroelectric materials.

Experimental Section
The calculations were carried out using the density functional theory (DFT) [58] as implemented in the VASP code. [59] The interaction of valence electrons with atomic cores is described by the projector-augmented wave (PAW) method, [60] as parameterized by the Perdew-Burke-Ernzerhof (PBE) functional. [61] It is demonstrated that, relative to other vdW correction schemes (such as many-body dispersion which is powerful in describing other systems [34,62] ), semiempirical DFT-D3 method [63] can more accurately describe the configuration of the ground sate of 1QL-In 2 Se 3 . [54] More calculation parameters and testing computational details are provided in Figures S1 and S6 in the Supporting Information. For 2D systems, the out-of-plane electric polarization is well-defined following the classical electrodynamics and calculated by the integration of local charge density times the coordinate in the out-of-plane axis over the whole supercell. [54]

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.